Construction of a Complete Genome Sketch of Bt S1028-1 Against Cylas formicarius

Zhang Jie ^1,2* , Wu Zhongqi ^1,2

1 Institute of Life Science, Jiyang College of Zhejiang A&F University, Zhuji, 311800

2 Cuixi Academy of Biotechnology, Zhuji, 311800

^* Corresponding author, jessi.j.zhang@foxmail.com

Abstract Bacillus thuringiensis strain S1028-1 showed high biological activity against Cylas formicarius. In this study, the de novo whole genome sequence of Bt S1028-1 was carried out by Illumina HiSeq 2000 and generated 1.64 Gb raw data. A total of 270 Scaffolds were obtained by quality control, filtering and SOAPdenovo short sequence assembly of the original reads data. The genome size was 5.87 Mb, N50 length was 356 721 bp, and the GC content was 34.99%. The genomic information of 17 Bt strains in NCBI database were selected as the reference, and the genomic sequences of Scaffolds were assembled and 3 replicators were obtained. Among them, the size of the pseudo-nuclear genome was 5.24 Mb and the GC content was 35.29%. The length of the genomic sequences of two plasmids were 517 kb and 116 kb, respectively, and the GC contents were about 32%. Based on the genome sequence and annotation information, the visualization sketch of Bt S1028-1 genome was constructed by GCviewer. The whole genome sequencing and construction of the draft of Bt S1028-1 provide the research basis for the functional study of this strain, and contribute to the distinction and identification of Bt S1028-1 insecticidal toxin protein.

Keywords Bacillus thuringiensis; Bt S1028-1; Whole genome sequencing; Genome assembling; Draft genome

苏云金芽孢杆菌(Bacillus thuringiensis, Bt)是一类包含许多变种的产晶体的芽孢杆菌，是革兰氏阳性细菌，普遍存在于大自然界中，如寄生在一些鳞翅目昆虫或者是植物表面等(Raymond et al., 2010)。Bt菌最初是日本人Ishiwata在1901年从病蚕中分离得到的，后来又有人不断的在其他地方发现此菌，1915年正式定名为苏云金芽孢杆菌。Bt在生长过程中会产生两大类毒素，分别是内毒素(伴孢晶体蛋白, insecticidal crystal proteins, ICPs)和外毒素，其中ICPs对昆虫具有一定的杀虫活性(Angus, 1956; Ibrahim et al., 2010)。经过研究发现，ICPs对鞘翅目、双翅目、鳞翅目和膜翅目等昆虫都有特异性的毒杀活性(de Maagd et al., 2001)，同时Bt是对环境友好的生防微生物，因此Bt是世界上应用最为广泛、用量最大、效果最好的微生物杀虫剂。

Bt S1028-1是由海南省热带农业资源研究所(HITAR)从广西大王岭自然保护区土壤中分离到的苏云金芽孢杆菌野生型菌株，该菌株对鞘翅目昆虫(甘薯小象甲(Cylas formicarius))和双翅目昆虫(Diptear)具有较高的杀虫活性(方宣钧, 2010)。甘薯小象甲属鞘翅目(Coleoptera)、象虫科(Curculionidae)昆虫，其成虫和幼虫均能危害甘薯(Dioscorea esculenta (Lour.) Burkill)，成虫咬食薯块、藤头、茎、叶柄、嫩梢、幼芽；幼虫在薯蔓、薯块内蛀食形成不定性的弯曲隧道，隧道内充满幼虫的排泄物，而且伤口易诱使致病菌侵入，使薯块发出恶臭和苦味，引发黑斑病、软腐病等，减低甘薯质量与产量。中国甘薯的年均种植面积为600万hm²，而甘薯小象甲危害造成的薯块产量一般损失为5%~20%，严重时高达30%以上，大多数疫区薯块为害率高达60%~100% (于海滨等, 2010)。Bt S1028-1作为对甘薯小象甲具有活性的生防微生物，具有其独特的优势，但是该菌株所含有的杀甘薯小象甲的毒素蛋白并不明确，为此，我们委托华大基因(中国深圳华大基因)对该菌株进行全基因组测序。

截止到2016年12月，在NCBI数据库中已收录了89个Bt菌株的基因组拼装和注释记录，包括36个全基因组数据、16个染色体数据、16个Scaffold数据和21个Contig数据(https://www.ncbi.nlm.nih.gov/genome/genomes/486?)，在这些已公布的Bt菌株基因组数据中，进行了全基因组gap closing得到全基因组数据的菌株相对较少，大部分的数据还是以Contigs或者Scaffolds的形式存在。随着全基因组测序技术的发展及测序成本的降低，越来越多的研究团队开展了Bt菌株全基因组测序的研究，通过深度测序获得质量更高的基因组数据，以期从基因组的水平对Bt个体和群体展开研究。本研究中，Bt S1028-1由华大基因进行Illumina HiSeq 2 000全基因组测序，然后对reads原始数据进行拼装得到了Scaffolds数据。并采用基于同一菌属类不同菌株间的序列相似性和局部最优搜索算法的生物信息学分析策略，从Bt S1028-1全基因组Scaffolds数据出发，构建了Bt S1028-1的基因组草图。这一研究结果将有助于开展Bt S1028-1菌株的基因组功能分析，对于挖掘新型Bt毒素蛋白具有重要的参考价值，为后续的毒素蛋白的功能分析提供了科学依据。

1结果与分析

1.1测序数据统计及Contigs/Scaffolds序列拼装

利用华大基因Illumina HiSeq 2000测序平台对Bt S1028-1菌株进行了de novo全基因测序，共得到1.66 Gb的原始数据，包含7 811 112个reads数据。通过对reads数据进行质量控制、过滤处理和SOAP de novo2短序列拼装处理，共得到270个Scaffolds，总长度为5 872 622 bp，N50值为66 528 bp，GC含量为34.99% (表1)。

表1 Bt S1028-1 Scaffolds的组装数据统计 Table 1 The statistic of assembled scaffolds of Bt S1028-1

1.2参考基因组的选择

目前NCBI数据库中具有基因组信息登记的Bt菌株共有89个，其中具有完整基因组信息的有36个菌株，只有染色体基因组信息的有16个菌株，还有37个菌株只具有Scaffold/Contigs信息。选取具有全基因组信息和只有染色体基因组信息的35个Bt菌株作为候选参考基因组，在ncbi-blast-2.2.27序列比对平台下，进行序列相似性分析。从序列一致性百分数、比对匹配的序列长度、序列比对错配数和空位比对上的碱基数四个方面，综合性分析基因组序列比对相似性结果，结果选择97-27-a、Al_Hakam-a、Bc601、BGSC_4AA1、BGSC_4C1、BMB171、Bt407、CT-43、CTC、HD-1、HD12、HD521、HD73、HS18-1、ST7、YBT-1518和YBT-1520a共17个菌株的基因组作为Bt S1028-1基因组测序拼接的参考基因组(表2)，其中菌株Bt407只有染色体基因组信息。利用DIYA基因组自动拼接工具中的拼接子程序diya-assemble_pseudocontig.pl，分别参照选择的17个已知Bt菌株基因组进行组装(图1)。

表2 参考菌株的基因组信息 Table 2 Genomic information of reference strains

图1 DIYA基因组拼装 注: 红色标记的为完成全基因组测序且血清型标记的菌株; 蓝色标记为完成全基因组测序但未血清型标记的菌株; 黄色标记为只有染色体基因组数据信息的菌株 Figure 1 DIYA genome assemble Note: Red labeled strain, complete genome sequencing and serological markers; Bule labeled strain, complete genome sequencing and not serological markers; Yellow labeled strain, only genome data information

1.3 Bt S1028-1基因组草图的构建

分别参照17个已知基因组进行组装，最终得到了1个拟核和2个质粒基因组，其中拟核基因组大小约为5.24 Mb，GC含量为35.29%；2个质粒的基因组大小分别为116 Kb、517 Kb，GC含量均在32%左右(表3)。利用Glimmer和GeneMarker软件对Bt S1028-1基因组分别进行开放阅读框(ORFs)和蛋白编码基因(CDSs)的预测，共得到6 169条ORF序列和5 916条CDS序列。利用GCviewer server以默认参数值对Bt S1028-1预测的基因组进行环形可视化草图构建(图2)。

表3 Bt S1028-1基因组信息 Table 3 Genomic information of Bt S1028-1

图2 Bt S1028-1预测基因组环形可视化草图构建 注: A: Bt S1028-1拟核基因组; B: Bt S1028-1预测的质粒1; C: Bt S1028-1预测的质粒2 Figure 2 Genome visualization graphical map of Bt S1028-1 Note: A: Nuclearoid genome of Bt S1028-1; B: Plasmid 1 predicted by Bt S1028-1; C: Plasmid 2 predicted by Bt S1028-1

2讨论

通常，全基因组测序原始数据中低质量数据的比例会在10%~20%之间。在本研究中，Bt S1028-1共得到1.66 Gb原始数据，其低质量数据约占原始数据的7.8%，这个过滤结果初步显示出Bt S1028-1基因组测序结果reads数据质量是比较高的，最终将92.2%的reads组装成了大小是5.86 Mb的Bt基因组。

Bt染色体为环状DNA，一般GC含量在32%~35%之间，基因组大小则是在2.4~5.7 Mb之间(Schnepf et al., 1998)。大多数的野生Bt内含丰富多样的质粒，而且质粒的DNA含量一般为细胞总DNA含量的10%~20%左右(Aronson, 1993)。本研究可知，Bt S1028-1包含一个染色体和两个质粒基因组，其基因组大小分别为5.24 Mb、517 Kb和116 Kb，GC含量分别为35.29%、32.55%和32.64%。从Bt基因组大小和GC含量来看，Bt S1028-1基因组的组装结果符合Bt菌株的基因组大小和GC含量的范围。通过GCview server对Bt S1028-1进行环形基因组的可视化发现，其基因组特征GC偏好性符合一般的Bt基因组的序列特征，这一结果说明Bt S1028-1基因组草图构建的合理性。

基因组组装结果的评价标准有N50大小、单碱基错误率、Contigo和Scaffold的错误组装等，而通常使用的评价标准同时也是组装结果评价的第一指标为N50，将Contig或Scaffold从长到短进行排序并累加，当累计和达到基因组序列50%是的Contig或Scaffold长度为N50。一般来说N50越长，表示组装得到的基因组序列可能更完整。当然N50长也不一定意味着组装结果一定可靠，因为当不相关的Reads或Contig错误的连接成Scaffold的时候也有可能达到较高的N50，因此在判断组装结果的时候要综合多方面考虑。本研究中，Bt S1028-1的N50为66 528 bp，低于细菌基因组N50参考标准300 Kb，初步判断Bt S1028-1的零碎序列较多，gap较多，但这并不能完全证明此次组装结果不理想，还要再进一步的对基因组组装结果进行验证。

在下一步的分析研究中，可以单独对Bt S1028-1所含有的各个质粒进行分离实验，进一步从分子水平检验Bt S1028-1基因组草图构建的正确性，或者进一步加大测序深度，对gap进行补缺失，对测序数据进行完善，以获得更为完美的结果。

Bt S1028-1的基因组测序和草图的构建，为Bt S1028-1菌株的功能性研究提供了研究基础，通过对Bt S1028-1的预测基因组进行编码蛋白预测和功能注释分析，有助于Bt S1028-1毒素蛋白的识别与鉴定，推动了Bt S1028-1菌株的应用开发进程。

3材料与方法

3.1 Bt S1028-1样本收集和基因组测序

Bt S1028-1是由海南省热带农业资源研究所(HITAR)从广西大王岭自然保护区土壤中分离获得的野生型菌株。由深圳华大基因研究院于2013年利用Illumina HISeq 2000测序平台完成了Bt S1028-1的全基因组测序。

3.2测序数据预处理

利用FastQC V0.6.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)对原始测序数据进行了质量评估，并对数据进行预处理。数据的过滤处理主要包括：a，去除reads质量值连续≤20的碱基数达到一定程度的readsmm (默认40%)；b，去除含N的碱基数目总和达到一定比例的reads (默认10%)；c，去除adapter污染(默认adapter序列与read序列有15 bp的overlap)；d，去除duplication污染。在后续的分析中，主要是针对clean data进行分析。

3.3 Contigs/Scaffolds的组装

运用SOAPdenovo2短序列组装软件(http://soap.genomics.org.cn/soapdenovo.html; r240)对clean reads数据进行组装，经多次调整后主要参数K设置为46，得到最优的组装结果。把clean reads比对到contigs上，根据clean reads的Paired-End和overlap关系，对组装结果进行局部组装和优化，最终形成Scaffolds拼装结果。

3.4参考基因组的选择

从公共数据库NCBI中下载已经完成全基因组测序工作具有完整菌株基因组信息的35个Bt菌株和具有染色体基因组信息的15个Bt菌株的基因组序列数据(Benson et al., 1990)，利用BLAST子程序makeblastdb分别对Chromosome序列和Plasmid序列构建比对数据库(Altschul et al., 1990)。利用ncbi-blast-2.2.27序列比对平台，进行各菌株基因组序列相似性分析，考察Bt S1028-1基因组与各菌株基因组的相似性，最终本研究选取了97-27、str. Al_Hakam等17个Bt菌株的基因组序列作为Bt S1028-1基因组测序拼接的参考依据。

3.5 Bt S1028-1基因组草图构建

利用有Paired-End关系的clean reads比对到Scaffold上，如果两个Scaffold之间存在潜在Paired-End关系，则推测这两个Scaffold可能存在前后关系，根据推测出的结果判断测序结果的可拼装性。并参考选取的NCBI数据库中已完成全基因组测序的17个Bt菌株的基因组序列，将组装得到的Scaffolds进一步分类组合形成Bt S1028-1菌株的pseudo基因组。再对Bt S1028-1基因组进行编码蛋白的预测和功能注释，利用CGView Server (http://stothard.afns.ualberta.ca/cgview_server/) (Grant and Stothard, 2008)完成Bt S1028-1预测基因组的可视化。

作者贡献

张洁是本研究的负责人，负责数据分析、初稿写作、修改和校对；吴仲琦参与数据整理。两位作者都阅读并同意最终的文本。

参考文献

Altschul S.F., Gish W., Miller W., Myers E.W., and Lipman D.J., 1990, Basic local alignment search tool, Journal of Molecular Biology, 215(3): 403-410.

https://doi.org/10.1016/S0022-2836(05)80360-2

PMid:2231712

Angus T.A., 1956, Association of toxicity with protein crystalline inclusions of Bacillus sotto Ishiwata, Can. J. Microbiol., 2: 122-131.

https://doi.org/10.1139/m56-017

PMid:13316606

Aronson A.I., 1993, The two faces of Bacillus thuringiensis: insecticidal proteins and post-exponential survival, Mol Microbial., 7(4): 489-496.

https://doi.org/10.1111/j.1365-2958.1993.tb01139.x

PMid:8384686

Benson D., Boguski M., Lipman D.J., and Ostell J., 1990, The national center for biotechnology information, Genomics, 6(2): 389-391.

https://doi.org/10.1016/0888-7543(90)90583-G

PMid:12134874

de Maagd R.A., Bravo A., and Crickmore N., 2001, How Bacillus thuringiensis has evolved specific toxins to colonize the insect world, Trends in Genetics Tig, 17(4): 193-199.

https://doi.org/10.1016/S0168-9525(01)02237-5

PMid:11275324

Fang X.J., Zhang W.F., and Xie L., 2010, Bt cry40a1 gene encoding insecticidal crystal protein that high toxic to Dipteran insects and their applications, China Patent, Application No.: 200910119992.9. (方宣钧, 张文飞, 谢柳, 2010, 对蚊子等双翅目昆虫有杀虫活性的Bt cry40Da1基因及其应用, 中国专利, 申请号200910119992.9.)

Grant J.R., and Stothard P., 2008, The CGView Server: a comparative genomics tool for circular genomes, Nucleic Acids Research, 36(Web Server issue): W181-184.

https://doi.org/10.1093/nar/gkn179

PMid:18411202 PMCid:PMC2447734

Ibrahim M.A., Griko N., Junker M., and Bulla L.A., 2010, Bacillus thuringiensis: a genomics and proteomics perspective, Bioengineered Bugs, 1(1): 31-50.

https://doi.org/10.4161/bbug.1.1.10519

PMid:21327125 PMCid:PMC3035146

Raymond B., Johnston P.R., Nielsen-LeRoux C., Lereclus D., and Crickmore N., 2010, Bacillus thuringiensis: An impotent pathogen? Trends in Microbiology, 18(5): 189-194.

https://doi.org/10.1016/j.tim.2010.02.006

PMid:20338765

Schnepf E., Crickmore N., van Rie J., Lereclus D., Baum J., Feitelson J., Zeigler D.R., and Dean D.H., 1998, Bacillus thuringiensis and its pesticidal crystal proteins, Microbiol. Mol. Biol. Rev., 62(3): 775-806.

https://doi.org/10.1128/MMBR.62.3.775-806.1998

PMid:9729609 PMCid:PMC98934

Yu H.B., Zheng Q., and Chen S.L., 2010, Biological characteristics and integrated control measures of Cylas formicarius, Hebei Nongye Kexue (Journal of Hebei Agricultural Sciences), 14(8): 32-35. (于海滨, 郑琴, 陈书龙, 2010, 甘薯小象甲的生物学特征与综合防治措施, 河北农业科学, 14(8): 32-35.)

Figure 1 New ICT based fertility management model in private dairy farm India as well as abroad